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ABSTRACT 



Two numerical methods of determining optimum 
paths in the problem of Bolza are presented. The basic 
theory of the methods is outlined and their essential 
characteristics are demonstrated by several specific 
applications to ship routing. 

The first method is characterized as one of differ- 
ential corrections employing methods developed by Professor 
Faulkner of the U. S. Naval Postgraduate School. The 
second method is termed a method of gradients employing 
calculation routines developed by Professor Bryson of 
Harvard. On the basis of the applications comparisons are 
made in the areas of ease of application and speed in 
producing results. 

I wish to express my appreciation for the guidance 
and encouragement of Professor Faulkner and for the 
programming of Mary Haynes of the Postgraduate School 
Computer Center. 
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INTRODUCTION 



The problem of obtaining maximum or minimum performance 
in a variety of controllable situations is undoubtedly as old 
as time. Usually the degree of success is determined by the 
efficiency with which some vital, limited quantity is ex- 
pended. 

Often the best solution cannot be defined, except as it 
develops by experience, even though the problem can be de- 
fined and stated mathematically. The magnitude of the 
analysis involved limits the feasibility of determining the 
solution. Modern high speed computing machinery has, for 
many problems, removed the latter limitation. As a result 
renewed interest has been shown and much progress made in 
the development and application of control theory to such 
problems of optimum control. 

This paper presents an outline of two practical nu- 
merical methods of solution of some general problems of this 
type. The class of problems is more specifically termed the 
problem of Bolza [l]. For purposes of identification in this 
paper the two methods will be called the 'differential 
method' and the 'method of the fundamental lemma' which will 
henceforth be abbreviated the D-method and the FL-method 
respectively. 

The D-method is based on the formulas for differentials 
due to Bliss [2], The numerical scheme for determining the con- 
stants of the solution by a Newton iteration is the work of 
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Faulkner [3], [4], His research papers and related studies 
are the source of both the theory and the numerical data 
presented for the D-method. 

In the literature the FL-method is also referred to as 
the 'gradient' and the 'steepest ascent' method. The calcu- 
lational procedure is similar to the proof of the fundamen- 
tal lemma and is apparently conceptually due to Courant [ 5 ]. 
The numerical scheme, which makes possible its application 
to the class of problems being considered, is the work of 
Bryson and Denham [6], A similar method is presented by 
Kelly [7], The theory and basic procedure for application of 
the FL-method, as reported in this paper, is that of Bryson 
and Denham. The contribution of the writer lies in the de- 
velopment of specific convergence schemes, left undefined by 
the authors [6], which make possible the programming and so- 
lution of two relatively simple ship routing problems thereby 
demonstrating several features of the method. Both problems 
have also been solved by the D-method and serve as a basis 
for limited comparison of the methods. 
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I 

INTRODUCTION TO PROBLEMS IN OPTIMUM CONTROL 

1.1 Formulation of a problem 

We desire the solution to the set of ordinary nonlinear 
differential equations 

( 1 ) 3C ^ ~ f L ( t ) y X 2 ( t) ^ 9 ♦ , , ^ X n (t) ) P}^ ^ ^ 1 9 9 } Pjjj ^ ^ ] 

i = 1,2, ... ,n; t Q t ^ T 

which is optimal, in the sense of maximizing the terminal 
value of a performance criterion 

g(x 1 , ,Xn,t) t=T , 

while simultaneously satisfying terminal constraints 

(2) h k (x 1? , x n»t)t=T =0 > k = l>2,...,Nj 0 < N < n. 

To effect these conditions we have at our disposal m con- 
trol or driving functions p (t), s = l,2,....,m. Their 

o 

definition for t Q 4 t ^ T is a major portion of the problem. 

* A dot ( ) over a variable is used to designate the deriv- 
ative with respect to t. 

If an integral is to be maximized we introduce an addi- 
tional state variable, x , , and an additional differen- 
tial equation 

*n+l = f n+l (x l’* ,, ’ x n’ t) 
where f n+1 I s the integrand of the integral. 
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It is assumed the functions f^ are defined over an 
open region and are class C". It is also assumed the func- 
tions g and h^ are independent and of class C* in the 
terminal region. The variables (x.^ ,x g , , . . . ,x n ) will be 
termed the state variables. All are functions of the in- 
dependent variable t. The terminal value, T, of t is 
either specified or otherwise determined from a terminal 
condition. The values (x-^ 

It will be convenient to represent these sets of varia- 
bles by matrices or by vectors. We define 



x 1 ( t) 




p^tf 




h^Xjt) 




f x (x,p,t) 


x 2 (t) 

• 


; p = 


p 2 (t) 

• 


; h = 


hg(x, t) 

• 


; f = 


f£^ X »P>^) 

© 


• 

*n (t) 




♦ 

p_(t) 

m 




© 

h N (x,t) 




© 



No distinction will generally be made between an n*l 
(column) or a lxn (row) matrix and a vector. If a quan- 
tity is best considered a vector, as in a scalar product, 
it will be designated, for example, pi'. 

Equations (1) and (2) then become 

(3) x = f(x,p,t) 

(4) h(x,t) t=T =0. 

We desire the solution of system (3) which maximizes 
the performance quantity g(x,t) subject to the con- 

w J* 

straints (4) . For purposes of introduction of terminology 
and some basic theory a brief review of the related elements 
of the Calculus of Variations is first presented. 

1.2 The Classical Calculus of Variations 

The problem formulated in section (1.1) is typical of 
those considered in the variational calculus - those of 



,x 2’** ,x n^to are assuine< ^ specified. 
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determining maxima or minima or, in general, stationary 
values (defined on page 8 ) of functionals such as g(x,t) 

by seeking the argument functions x(t) and p(t) for 
which the functional assumes the stationary value or ex- 
tremum in question. 

An outline of seme phases of the classical method of 
the variational calculus, as applied to a greatly simpli- 
fied optimum problem, will serve to introduce terminology 
and to present essential concepts which are the basis for 
both methods of solution being considered. The method may 
be termed 'indirect' in that it is based on the reduction 
of the variational problem to differential equations the so 
lution of which, subject to the boundary conditions, is the 
solution to the problem. The two methods considered in 
this study may be termed 'direct' in that they consist of 
construction of a sequence of functions that converges to 
the desired function. 

Following Courant and Hilbert [S], we consider the 
simplest problem of the variational calculus. Find the 
function y(t) which is such that y(t Q ) = A and 
y(t^) = B, Fig. 1, and the functional 



is maximized. Any function y = y(t) which meets the end 
conditions, is of class C, and has a piecewise continuous 




T 



o 



• 5 



t 




hereafter we will' refer only 
to maxima) the functional 
F[y( t) ]. Suppose we have 



Fig. 1 An admissible 
variation of an admissible 
curve 



found the desired extremal function y = y(t). Let us con- 
sider also an unspecified new function n(t) which is de- 
fined for t 4 t < t, , possesses the necessary derivatives 
and continuity, and is such that n(t 0 ) = h(t^) = C, but is 
otherwise arbitrary. The function 



y = y(t) + c t) ( t ) = y(t) + 6y(t), 



where e is a parameter, represents a one-parameter family 
of admissible functions which contains the solution y(t) 
when €=0. The quantity 6y(t) = e*n ( t) is the vari a tion 
of the function y = y(t). For sufficiently small magni- 
tudes of c the varied functions y(t) lie in an arbi- 
trarily small neighborhood of the extremal, y = y(t), and 
the integral 

F(y) = F[y( t) + eri(t)] 
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ay be regarded as a function i(c) which t ust have a h axi- 
mum at e = 0 and therefore requires 'I'(C) = 0. 

Kow 



implies 



y = y( t) + en( t) 



y'(t) = y * ( t) 4- en'(t) 



and we may write 



$(e) = / f[t, y + € rt, y'+en'] dt 

*o 

For a maximum it ir necessary that 



(5) 



§ * Co) 



di , 

■ L f y * 3tV ] dt =° 



for all n( t) meeting the conditions previously imposed. 
The fundamental lemma of the calculus of variations, 
[8] page 185, applied to (5) requires 



( 6 ) 



•p _ ;L.f = r 

"y dt y 1 c ’ 



or equivalently, 



(7) 






x y'y 



+ V 



* f 

y'y' 



+ f 






y ' x 



- f y = c 



This differential equation was first discovered by 
Euler in 1744, [9j page 22 , and is commonly referred to as 
the Euler equation. A solution to the Euler equation is 
called an extremal. 



Its validity is a necessary condition 



for the existance of an extremum^S] page 185. Being a 
differential equation of the second order, two arbitrary 
constants must be determined for a solution which meets the 
boundary conditions, i.e., to make an extremal an admissible 
extremal. Every solution of Euler's equation is an ex- 
tremal of the maximum problem, [8]. Courant describes (6) 
the variational derivative of F with respect to t and 
terms its role here as analogous to that of the gradient in 
ordinary maximum problems. 

In the n+1 dimensional case, with coordinates 



t, y x (t), y 2 (t), . . . , y n (t) 



and the functional 

_t, 

( 8 ) 

■ A o 

the Euler equations take the form 



t l 

1 ~ ft f(t ’ 7 2 ’ **’’ 7n ’ y 2 ,# 



•> y n ) 



dt f y n ' ” = " f v_" " f v ” 



dt y 



‘dt y ' 



1 J 2 J 2 "n "n 

a system of the second order equal in number to the number 
of unknown functions y^(t) to be defined. Any n+1 space 



C VC 



y l = y l (t) > y 2 = y 2 (t) ’ » y n = y n (t) 

v/hj ch is a solution to the system (9) is an extremal and 
furnishes a stationary value of the functional (8). 
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The necessary condj tion of Euler is usually called the 
first necessary condition. There are also necessary con- 
ditions of Veierstrass, Legendre and Jacobi as presented, for 
example, in [l] chapter 1. These are not discussed here; 
the fundamental problem confronting an applied mathematician 
is that of finding a likely candidate for a solution. 

1.3 Introduction to direct numerical methods 

The problem as stated in section (1.1) is similar to 
that resulting in system (9). The state variables x(t) of 
our functional g(x,t) are defined by a set of n first-order 
differential equations and one or more control functions, or 
parameters, p(t); t is the independent variable. Solu- 
tions to the related Euler equations will, to the extent of 
meeting the first necessary condition, define the functions 
p(t) such that the resulting path x(t) = C will impart 
at least a stationary value to the functional g(x,t)^,. 

The proper choice of the arbitrary constants of these so- 
lutions will result in admissible extremals if they are 
chosen such that the terminal conditions h(x,t)^, = 0 will 
also be met. It will be assumed that the initial conditions 

x(t). . are specified constants. 

r o 

Various indirect approaches to such a solution exist. 

A bibliography to indirect method studies is presented in 
the survey paper [10]. Analytical solutions of the Euler 
equations involved usually require idealizing assumptions 
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which limit their applicability in practical situations. 
Under more realistic assumptions a numerical attack on 
these equations is required. Some direct methods are pre- 
sented in [8] page 174. In the following sections the 
essentials of two additional practical methods are out- 
lined. As with all methods of solving Bolza's problem, 
Lagrange multipliers are employed producing equations 
equivalent to the Euler equations in terms of solutions 
to an adjoint system of differential equations. The con- 
cept was developed by Bliss [2] in connection with 
differential corrections in ballistics. 

The D-method uses only extremals. Then, by an iter- 
ative correction procedure, obtains an admissible extremal 
out of the family of extremals. 

The FL-method produces an extremal out of any likely 
solution to the original equations; part of the problem is 
in a subroutine to approach an admissible solution. The 
solutions to the Euler equations, producing the extremal, 
are approached iteratively by a gradient, or steepest- 
ascent technique. 

Since all methods of solution involve solutions to 
the adjoint equations this concept will be outlined be- 
fore the details of the two methods are presented seper- 
ately . 

1.4 The adjoint equations 

The problem as stated in section d.l) was to determine 
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the solution to the system 



(3) 



x = f(x,p,t) 



such that the terminal constraints 



(4) 



h(x,t) T = 0 



are satisfied and the functional g(x,t)^ is maximized. 

Of primary interest is the relationship between vari- 
ations fip(t) in the control variables p(t), at any 
value t, and the resulting terminal variations 6x,j of 
the state variables. We also need the relationships 
between 6p(t), fih^ and fig T . 

Looking at (3) in component formed), the variational 
relationships for the state variables can be expressed as 





i 1* 2. «•«.. n . 



In matrix notation this is 



(10) fix = F fix + G fip 
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where 



a xt 



fl£i 

* 9x £ 



m, 

0Pi 



Mi 

a pz 

m 



F(t) = 



; G(t) = 






af 

5? 



St . . 56. 



nxn 



8P 



8 Pi 



m 



J nxm 



6p = 



6p x (t) 



fix. 



fix = 



ep m (t) 



8i n 



We will use the prime symbol to denote the transpose of a 

matrix, e.g., fix' = (6x 1 , fix 2 , fix^. 

Following Bliss [2] the system of linear equations 
adjoint to (3) is defined to be 



(12) A= -F'(t)A 



where 

A'(t) = [x^t), x 2 (t).. x n (t)] 



is a matrix of n Lagrange multipliers, and 



A'= t i v 




o © © e o © y 




]. 



The solutions of systems (3) and (12) are related by 
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citLA'fix] = AG 6p; 



hence * 







Since we are assuming x t fixed fix^ = 0 in (13) and we have 



This is the differential formula relating the terminal change 
in the state variables and the variations in the control 
variables. It is termed the Fundamental Formula by Bliss [2] 
It is also known as the Fundamental Adjoint Formula and as 
the one-dimensional form of Green's Theorem [ll], [12]. 

1.5 Applications of the Fundamental Adjoint Formula 

Some indication of the usefulness of the fundamental 
formula is displayed in the following observations. 

1) If A is chosen as the specific solution to the 
adjoint system (12) such that at t = T 



o 



o 




T 



(14) 



o 



7 \ t = V^d) = (1,0,0,...,0), 



V^x-^ refers to the gradient of x-^ in x space. 



We assume no corners, [9] page 8, exist 
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let us call this solution /^ ( t) , (14) then becomes 

6x,(T) = s/a: G 6p dt , 

*0 

yielding the terminal variation of x^ due to variations 
6p(t). Similarly choosing A i (t) such that 

A.(T) = V x . ( T ) = (0,0,. ..,1, ,0) 

-L 

i = 1,2, ,n 



gives T 

(15) 6 x . ( T ) = f. G 6p dt . 

1 z o 1 

2) If A is chosen as the solution to the adjoint 

system such that A (T) = Vg(T), (14) implies 

6 * 

r T , 

(16) 6g(T) = (Vg * 6x) T = J t /[' G 6p dt 

* 1 ^o g 

giving the first variation of the performance function in 
terms of the parameter variations for t Q ^ t T. 



3) For the constraint functions we likewise have 



(17) 



T 

6h.(T) = J' f\* G 6p dt j j = 1) •••••> k* 



‘o 3 



4) If we were to ignore the problem of meeting con- 
straints and look for a solution which gives a stationary 
value for g(T), we would require the vanishing of the first 
variation, ( see appendix of [4]) 



(18) 6g(T) = / /[* G 6p dt 



v g 
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That is 



which requires A^ G = C. 



Ui >A 



2 ’ 



»A n ) g 



8p l 



ep 



■n 



~1 

8p n 



& 

dp: 



m 



= C 



which represents the m equations 



(19) 



K • S5. ■ o. i ■ i, 



,m. 



‘g QPi 

Equations (19) with the adjoint differential equations 
(12) constitute what is generally termed the Euler-Lagrange 
equations [4], [13], 

In equations (18) the matrix /\‘ G is clearly an in- 

o 

fluence function giving for any time t the effect that a 
unit impulse fp } at that time 5 would have on the perform- 
ance function g(T). In that regard if we write /^G as 

o 

a vector, say A(t), (18) then can be written 



» - £ 



*g(T) = L A(t) • «p dt 
p o 

and for maximum change in the performance quantity g(T) 
we would choose Cp parallel to A(t) for t ^ t ^ T. 

It should be noted that A(t) and fp(t) are m-dimensional 
vectors in parameter space (p). Courant, [8] page 223, re- 
fers to vector A(t) as the gradient of the functional 
g(t) in function space. When the optimum solution or path 
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in function space is obtained this gradient will vanish. 

In the case of our general problem we have, in addition 
to the requirement of maximizing the functional g(x,t)^, the 
necessity for meeting terminal constraints h(x,t) T = 0. We 
will choose the vector p(t) in parameter space such that 
the path meeting the constraints yields the largest possible 
value for g(T). We now consider two methods of solving the 
problem. 

We have called g(x,t)^, a functional. It is a function 
which becomes a functional since (x,t) represents the endpoint 
of a curve^ satisfying the given differential equations with 
initial point x(t ) given. The function g has become a 
functional in the sense that the point (x,t) depends upon 
the curve. 
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II 



TV/O NUMERICAL METHODS 



2.1 The differential method 

Recalling the fundamental differential formula 

(14) [A6x ] T ~ J G 6p dt , 

o 



the D-method is based on the Euler equation which states 
that if the path is to furnish a maximum to some functional 
g(x,t)rp it is necessary that the coefficient Ag of 6p 
in (14) vanish for some solution A of the adjoint system 
for all t t T. This requirement, Ag = 0, expanded 
in terms of (11), implies that for all t Q <C t ^ T there 
exists a A(t) and x(t) such that 



( 20 ) 



A(t) 

A<u 




= 0 





*3f _ 

ap 



The problem becomes: determine A(t) and hence the control 
functions p(t) such that (20), (3), (12) and (4) are satisfied. 

Assuming first that T is specified we proceed to out- 
line the method of determining the A matrix. On page 14 
a set of n solutions of the adjoint system (12) was 
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. 























developed such that 



Aid) -VjEjd), i = 1, 



,n. 



The linear combination 



(21) A(t) - c^ + C 2^2 + 




of these n solutions is also a solution to (12). 

As a first approximation to the admissible extremal 
solution desired we estimate a set of constants c^ for 
(21) and determine the associated extremal by employing 



an admissible extremal as it will not meet the terminal 
constraints (4) . 

On the basis of the error in meeting the terminal con- 
ditionsy corrections dc^ for the constants in (21) are 



with (3) to obtain a new solution, will, if the method con- 
verges, give a more nearly admissible extremal. The process 
is repeated until some convergence criterion is satisfied. 

If the terminal T is not specified it too will be es- 
timated for the first solution and a correction dT will be 
produced with the dc^ at the end of each iteration. 

A more detailed analysis of the basis and the calcu- 
lation procedure is presented in appendix (A). 



this A(t) in system (20). This will in general not be 



determined. 



The resulting improved A(t), when employed 
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2.2 The method of the fundamental lemma 

This method is concerned directly with the construction 
of the set of m control variables p(t), for t Q ^ t ^ T, 
which produce the solution to (3) that satisfies the terminal 
conditions (4) and maximizes the performance quantity 
g(x,t) T . Essentia ll y the procedure is "to produce any reason- 
able trial solution by the choice of a nominal control pro- 
gram p(t). The resulting solution will in general be 
neither admissible nor an extremal. That is, the con- 
straints (4) are not satisfied and the Euler equations 
(19), a necessary condition, are not satisfied. A better 
solution can therefore be constructed. A variation 6p(t) 
is then produced such that the Euler equations are more 
nearly satisfied and the solution's deviation from that 
of an admissible extremal is reduced. This variation ap- 
plied to the nominal p( t) produces a better control pro- 
gram ptt) = ptt) + 6p(t). The procedure might be de- 
scribed as iteratively stepping a control vector p(t), 
defined in parameter space for t Q < t ^ T, in the direc- 
tion of steepest ascent toward meeting constraints and 
toward meeting the first necessary conditions of an ex- 
tremal. This direction is that of -V h and V g in 

P P 

parameter space, where we use the subscript (p) to denote 
parameter space. As an admissible solution approaches an 
extremal V^g will tend to zero [6]. The direction of 
6p, for each point of a nominal solution, will be determined 
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by the direction cf the gradients at that point. The 
magnitude of 6p at each point will be best chosen pro- 
portional to the magnitude of the gradients at that point 
in parameter space. 

In appendix B the desired parameter correction program 
p(t) is derived and defined for t 0 ^ t < T as 
( 22 ) 

6 P (t) =±o'(/i ff A. i;i I*.> 



(ds) 2 - dh'i:?; dh , 

gn' M hn A hh ■ L hg '/“ 7 ~ 7* + G ^ hn Ihh dh * 



gg 



“hg hh hg 



The plus sign is used if g(x,t)^, is being maximized. 

One phase of the calculation procedure as outlined by 

Bryson [6], and summarized in appendix B, was found to be 

somewhat indefinite. It is that of determining the arbi- 

o 

trary constant (ds) which governs the magnitude of the 
step size §p(t) by limiting the integral of these mag- 
nitudes for t Q < t < T by the relationship (B-7) . It 
seems apparent that the larger the magnitude (ds) , within 
the limit which allows the linearization (10), the more 

rapidly our solution process will converge to the admissible 

2 

extremal. Accordingly the first restraint on (ds) will 

2 

be imposing an upper limit (ds o ) as estimated from the 
nature of the problem in conformance with (10). 

Further it is obvious that repeated application of 
large magnitude variations ^p(t) will not allow fine ad- 
justments toward the admissible extremal, as the solution 
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o 

is approached, unless the magnitude (ds) is further ma- 
nipulated, We therefore desire an automatic scheme which 
will rapidly step the nominal solution toward the ultimate 
solution but will adjust itself to finer increments in the 
vicinity of the admissible extremal. The following dis- 
cussion of Bryson's procedure is to develop a basis for the 
design of such an automatic convergence scheme. 

The expression (22) defining the vector 6p consists 
of two components which we define 6p^ and 6p » We 



g 



designate 



(23) 



6 fh = G ' A ha I hh dh 



since it is essentially a vector in p space whose magni- 
tude is proportional to the magnitude of the error in meet- 
ing the constraints resulting from the nominal solution. 

has the direction, essentially, of -Vph. Accordingly 
6p"k is a component of 6P which produces a more nearly 
admissible solution. The desired correction toward an 
admissible solution may be large and may require larger 
values of Up I than the upper limit, previously imposed, 
will allow. We give construction of an admissible so- 
lution first priority, however, and up to the limit (ds o ) 
demand that gp be designed to meet the constraints. 

The first term of (22), involving the radical, is 
essentially a vector in p-space whose direction is normal 
to V h and, as nearly as meeting constraints will allow, is 
parallel to V g. It specifies a component 6p which causes 

r* o 
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progression toward maximization of g(x,t)^, without dis- 
turbing the progress toward an admissible solution made by 
6p^. The term 

/(ds) 2 - dh' I{£ dh 

2 

is what remains of magnitude (ds) for use as magnitude for 

«. 2 

a vector normal to op^. If 6p h demands all of (ds) 

this radical is zero and 6p = 6p^. Progress will be made 
toward an admissible solution and the functional g(x,t)^ 

will increase or decrease according to whether or not Tp^ 

has a component in or opposite to the direction of Vpg. 

Accordingly, the scheme for the choice of the magni- 
' 2 

tude (ds) will be to first provide,' within the upper 
limit, a magnitude 

I °’A hn C I 

V 

for attaining an admissible solution, and second, provide 

for a component toward the extremal, adjusted according to 

the estimated deviation from an extremal. We define this 

latter magnitude (A), the magnitude of 6p . We then re- 

& 

quire 

(24) (ds 0 ) 2 (ds) 2 = a 2 + |G’ A hn dh | 

The initial choice and automatic adjustment of the 
magnitude (A) is another problem. We want it as large as 
possible initially. This will be controlled by (24). 
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We want smaller magnitudes as the ultimate solution is 



approached, automatically adjusted to the need. The con- 
vergence criterion suggested by Bryson [6] v/ill now be em- 
ployed. He develops the fact that on an admissible solu- 
tion the denominator of the radical of (22) is the direc- 
tional derivative or gradient 



in p-space. This must decrease as the extremal is approached. 
This quantity is, however, in general greater than zero. It 
is reasoned that if a relatively large (A) is repeatedly 
applied to an admissible solution it will rapidly decrease 
the quantity dg/ds of (25) and then will ultimately cause it tc 
pass through a minimum, greater than or equal to zero. If 
at this point the magnitude (A) is not reduced subsequent 
solutions will be over-corrected toward an extremal; trajec- 
tories which oscillate about the solution to the problem. 

We use this characteristic to trigger a reduction of (A) 

* 

to some fraction, say 0.2 A, and step from there toward 
the extremal in smaller steps. The minimum value for dg/ds 
attained this time will be smaller and the same behavior will 
signal the need for further reduction in (A). The process is 
repeated until some lower limit on (A) or dg/ds, a conver- 
gence criterion, is attained. 

A flow chart of a convergence scheme is presented as 
appendix C. 



( 25 ) 
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APPLICATIONS AND COMPARISON 



3.1 A basis for comparison 

The performance of both methods was investigated for 
two ship renting problems. Of primary interest were: 



(a) the problem of choice of starting para- 
meters that result in convergence to a 
solution. 

(b) the domain of these parameters that pro- 
duces a solution. 

(c) the range of problem variations that can 
be solved with an acceptable set of these 
parameters. 

(d) the rate of convergence to a solution. 



For each application the performance criterion was 
defined as the functional 



which could be interpreted, for example, as a probability 
of detection or as a measure of accumulated radioactive 
fallout in an area where the intensity is defined by the 
time and position function g(x,t). This is essentially 
the problem considered by Faulkner [4] and Cook [14]. 
Their programs for the D-method were employed as a source 



(26) 




(x - t/100) 2 + 4v 2 
2 



] dt 
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of comparison data for that method. 

More specifically the problems considered are: 

Problem I: Determine the route from a point within an 

area where g(x,t) is defined, to any point 
at which g(x,t) is a prescribed tolerable 
maximum, such that the accumulated fallout 
g(x,t) T is minimized. 

Problem II: Determine the route between two prescribed 
points in the area where g(x,t) is defined 
such that g(x,t) T is minimized. 

For the D-method the calculation procedure of page 18 
for these problems is reduced to the specification of a 
pair of constants defined by Cook [14] as (X Q ,u 0 ) and are 
the two components of the initial adjoint vector, which are 
unknown at the outset of computation. 

For the FL-method the starting parameters are, as for 
any problem, a nominal p*(t) and a magnitude (ds o ) . For 
these problems p*(t) is a bearing angle or heading and 
could be any constant direction p*(t) = k radians, rough- 
ly away from the point of maximum g(x,t), 'ground zero'. 

From (B-7) an acceptable value for (ds Q ) is determined by 

2 2 

(ds^) = (6 p) T, where T is approximately the total time 

required and 6p is the change in bearing, at any point, 
that would reasonably meet the linearity requirements of 
(10). As an example; estimating T as 100 with an allow- 
able 6 p of 0.5 radians would give: 
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(ds o ) 2 = (0 . 5) 2 ( IOC) = (b) 2 . 

3.2 Problem I 

We consider the specific case of problem I with start- 
ing point (x ,y ) specified as (0.3, 0.1). For the FI..- 

methcd the choice dS^ =5 is reasonable and a logical 

o 

choice for p*(t) is p*(t) = v/2 for 0 < t ^ T or very 
roughly 'away' from g(max) at (0,0). For the D-rnethod a 
reasonable choice (with the advantage of hind-sight) for 
(X o , u 0 ) is (0,-100) whj ch is a vector roughly in the direc- 
tion of the negative gradient of g(x,t) with a magnitude 
approximating the final time T« 

In an effort to experimentally determine the intervals 
of convergence an interval of values containing these esti- 
mated starting parameters was explored for both methods. 

The results are summarised in Table I and Figures 2 and 3. 

For the FL-method the process converged for 
0.1 < p*(t) ^ 3.2 radians. The starting parameter is 
therefore any bearing angle roughly 'away' from (0,0) within 
about a 180° interval. In Fig. 2 path A depicts the op- 
timum route. For values in the vicinity of 3v/2 the local 
minimum at the opposite side of the cloud (path B in the 
sketch) is developed. If headed in the direction p v (t)= v y 
however, the local maximum will not be produced but the 
nearest minimum (path A) will be generated, while p#(t)= 3.6 
results in convergence to pathB. 
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Table I 



Starting parameters for problem I. 



D-Method 


FL-Method 


>0 


MO 


Convergence 


p*(t) 


Convergence 


90 


-10 


yes 


0.1 


yes 


80 


-20 


yes 


0.2 


yes 


60 


-40 


yes 


0.4 


yes 


40 


-60 


yes 


0.8 


yes 


20 


-80 


yes 


1.2 


yes 


0 


-100 


yes 


1.6 


yes 


-20 


-80 


no 


2.0 


yes 


-40 


-60 


no 


2.4 


yes 


-60 


-40 


no 


2.6 


yes 


-80 


-20 


no 


3.0 


yes 


-90 


-10 


no 


3.2 


yes 


0 


-50 


no 


3.6 


(path B) 


0 


-200 


yes 






60 


40 


(path B) 






80 


20 


m 






0 


100 


it 






-60 


40 


no 






-80 


20 


no 
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V 




Figure 2 FL Convergence Interval for p*(t)=k 



The number of iterations required for a 'solution' de- 
pends on the accuracy demanded. For values within the inter- 
val of convergence the poorest choice for p*(t) requires 
about ten percent more iterations than required of the 
best constant choice. No great penalty is associated with 
rough estimation of the starting parameter, for a close 
approximation to an admissible route is produced rapidly. 
Computer time is of the order of 45 seconds for 40 itera- 
tions. The first 15 iterations produced a very accurate 
solution. 

For the D-method starting parameters were explored in 
the intervals -100^ X c ^ 100, -100 < u 0 100. For pur- 

poses of probing the influence of magnitude (X 0 ,u o ) =(0,-50) 
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Table I sum- 



ancl ((' , -f('O) were also tried. Figure 3 or/ 
marine the results. In each case for which a solution 
existed convergence was relatively rapid requiring about 
10 iterations and computer time in the order of 10 seconds 




for the particular convergence criterion employed. 

Continuing with problem I we explore the domain of 
starting points ( x 0 >y 0 ) for which a given set of starting 
parameters produces a solution. For practical purposes it 
would be desirable that a given set of parameters func- 
tion for a large domain of starting points under the radio- 
active cloud. Accordingly solutions were attempted for a 
number of starting points in an interval bounded roughly by 
the left and right extremeties of the intolerable fallout 
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area. The input parameters were maintained at (X o ,u o )=(0,-100) 
for the D-method, and p*(t)=1.0 for the FL-method. 

Figure 4 and table II portray the results for eight start- 
ing points. Both methods produced solutions for each run. In 
figure 4 the terminal points are associated with the corre- 
sponding initial point by a dashed line simulating the route. 

As table II indicates, the coordinates of the terminal points 
agree for the two methods to the third or fourth significant 
figure. Values for final time (T) and the performance g like- 
wise agreed to at least four figures. Extended results for them 
are given only for run D , 

In regards to the FL-method, table III is a printout of 
the terminal values of the iterations toward a solution for 
run D of figure 4 and table II. It will serve to display 
several features of the FL-method as applied to a relatively 
extreme choice of the starting parameter p*(t). Specifically 
the starting point is (0.3, 0.1) and p*(t) was, for probing 
purposes, chosen as p*(t)=0.1. This is an obviously unwise 
choice in as much as the cloud and the ship have comparable 
speeds and, for this case, are almost parallel. Consequently 
the terminal time (T) for the nominal run is about ten 
times the value on the extremal. A nominal solution is pro- 
duced, however, and from it relatively rapid convergence to 
an extremal is portrayed. 

The parameter (dS Q ) was chosen as (6.0). The choice of 
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Fig. 4 Allowable Starting Points for (> 0 ,u 0 ) and p*(t) fixed 



Table II Allowable Starting Points (x ,y ) for Problem I 
(> 0 ,u 0 ) = (0,-100); p*(t) = 1.0 



Start point 
(x o> Jr o ) 


X(T) 


Y(T) 


FL 


D 


j FL 


D 


A (1.2, .1) 


1.4092 

i 


1.4103 


| 

| 1.0511 

j 


1.0510 


B (1.0, .1) 


1.1379 

l 


1.1388 


| 1.0702 

1 


1.0702 


C (0.7, 0.1) 


! 0.7266 

I 


0.7223 


j 1.0661 


1.0662 


D (0.3, 0.1) 


0 . 1775 


0.1778 


1.0074 


1.0075 


E (-0.7, 0.1) 


-1.1066 


-1.1071 


0.6128 


0.6125 


F (-1.0, 0.1) 


-1.4175 


-1.4179 


0.4406 


0.4401 


G (-1.2, 0.1) 


-1.5879 


-1.5882 


1 

0.3325 

; j 


0.3320 


H (-1.4,0.!) 


-1.7292 


-1.7293 


0 . 2454 | 


0 . 2452 



Typical values for final time T and performance g(T) 

for point D: Tp L = 91.6399, T^ = 91.6417 

g( T) p j -= 49 .0500 , G(T) d = 49.0501 
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Table III Problem I, FL-Method, Run D 



X o = 0 . 3000 


' Y o = 


0.1000 


S = 6.00 


p X (t) = 


0.10 


T 


X(T) 


Y(T) 


g(x,t) T 


D 


ds 


967.1896 


9.9236 


1.0656 


487.8577 


17661.3315 


6,0000 


256.7626 


2.6762 


1.0716 


127.6031 


2228.5757 


6.0000 


124.8307 


1.0787 


1.0696 


62.3428 


487.7631 


6.0000 


95.2305 


.3507 


1.0300 


49.8881 


79 . 2727 


6.0000 


98.7530 


.1648 


.9910 


53.9286 


275.8876 


1.8000 


94.0974 


.1323 


.9939 


50.3492 


132.3472 


1.8000 


91.8374 


.1091 


.9938 


49.1760 


29.3284 


1.8000 


92.2926 


.2244 


1.0146 


49.3411 


64.2659 


.5400 


91.8649 


.2062 


1.0121 


49.1077 


23.8514 


.5400 


91.6341 


.1808 


1.0080 


49.0658 


15.8866 


.5400 


91.7693 


.1826 


1.0081 


49.1008 


28.9107 


.1620 


91.6957 


.1801 


1.0078 


49.0645 


15.1434 


.1620 


91.6437 


.1775 


1.0074 


49,0503 


1.5992 


.1620 


91.6199 


.1759 


1.0072 


49.0569 


11.0875 


.0486 


91.6262 


.1765 


1.0073 


49.0528 


7.1449 


.0486 


91.6340 


.1771 


1.0074 


49.0505 


3.2402 


.0486 


91.6426 


.1776 


1.0074 


49.0501 


.7613 


.0486 


91.6293 


.1766 


1.0073 


49.0503 


3.1283 


.0146 


91.6331 


.1769 


1.0073 


49.0501 


1.9116 


.0146 


91.6372 


.1773 


1.0074 


49.0500 


.7058 


.0146 


91.6416 


.1776 


1.0074 


49.0501 


.4999 


.0146 


91.6376 


.1773 


1.0074 


49.0500 


.7230 


.0044 


91.6388 


.1774 


1.0074 


49.0500 


.3542 


.0044 


91.6399 


.1775 


1.0074 


49.0500 


.0208 


.0044 


91.6395 


.1774 


1.0074 


49.0500 


.3385 


.0013 


91.6396 


.1774 


1.0074 


49.0500 


.2254 


.0013 


91.6397 


.1774 


1.0074 


49.0500 


.1132 


.0013 


91.6399 


.1774 


1.0074 


49.0500 


.0149 


.0013 


91.6398 


.1775 


1.0074 


49.0500 


.0327 


.0004 


91.6399 


.1775 


1,0074 


49.0500 


.0116 


.0004 


91.6399 • 


.1775 


1.0074 


49.0500 


.0062 


.0004 


91.6398 


.1775 


1.0074 


49.0500 


.0189 


.0001 


91.6399 


.1775 


1.0074 


49.0500 


.0106 


.0001 


91.6399 


.1775 


1.0074 


49.0500 


.0031 


.0001 


91 . 6399 


.1775 


1.0074 


49.0500 


.0040 


.0000 


91.6399 


,1775 


1,0074 


49.0500 


.0000 


,0000 



Computer times 44 seconds 
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any number 0 < ds o < 20 would also produce solutions. Actually 
values nearer 20 would hasten the convergence while values 
smaller than six would require a greater number of iterations. 

Table III also portrays the operation of the convergence 
scheme discussed on page 23. Since in this case there are no 
constraints , A in (24) is (ds). As table III indicates, dS 
stays at its initially chosen value of (6.0) for four itera- 
tions. By this time the extremal has been closely approxi- 
mated not withstanding the difficult nominal route chosen. For 
the fifth iteration dg/ds has crossed through its minimum and 
the reduction of (ds) to 0.3(dS) has been switched. The pro- 
cess then continued until an extremal was crossed again, caus- 
ing a further reduction in step size. In this case the pro- 
cess was halted when the condition dg/ds < 10'^was attained. 

3.3 Problem II 

The two methods were less directly compared for problem II. 
Experimental probing was directed primarily at weaknesses of 
the convergence scheme devised for the FL-method. The first 
specific problem of type II was: Determine the optimum route 

from (x o ,y o )=(0.03,0.1) to (x f ,y f )=(l. 0,1.0) with c = 0 
in (26). An admissible solution, simulated by route (A) in 
figure 6, was closely approximated in the third iteration. 
Convergence to an admissible extremal, however, was not re- 
alized. Longer and longer admissible routes, (B) , were devel- 
oped, each better than the previous from the standpoint of 
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decreasing the performance criterion, g . It then became appar- 
ent that the solution to this problem does not exist. The quan- 
tity g clearly would be minimized by escaping from under the 
radioactive cloud to some very distant point, returning after 
the cloud has moved away. The FL-method developed a sequence 
of solutions which approached this solution. The D-method for 
this case returned an' error* printout in as much as one of the 
input parameters is a relatively close estimate of the final 
time T. This value being undefined the process stopped due 
to too large an error in the estimate of the magnitude of this 
parameter. 

The FL-method was then applied to the same problem with 
c = 0.1 in (26). An admissible extremal resulted in a reason- 
able number of iterations. 

The importance of the careful choice of the stopping con- 
dition fj(x,t)^,= 0, for the FL-method, was demonstrated by 
problem II for the route (0.03,0.1) to (0.0, 1.0) , s .figure 6, 

y ' B ^ 

/ \ 




x 





< i / 

\ If 



i(t)=k 



X 

V 




X 



/ 



Figure 5 Problem Ila 



Figure 6 Problem Ilb,llc 
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with n(x,t) T = 0 chosen as (x - Xp)^, = 0. The resulting 
situation is that the stopping constraint is nearly parallel 
to the route making T, and thus x(T) and g(x,t)^ particularly 
sensitive to variations in the control variable. Consequently 
difficulty arose in that the constraint (x - x p ) T = 0 was 
difficult to meet with the accuracy desired. The choice of the 
line (y - y F ) T = 0, or the circle through (x F ,y F ) with 
center (x 0 ,y 0 ), would alleviate the situation. In general 
n(x,t) T = 0 and g(x,t) T = 0 should be as nearly parallel 
as possible. 

The difficulty arising from too large a value for A Q was 
demonstrated for the FL-method by the solution of problem II 
for the route (0.03,0.1) to (-1,0.1), simulated by curve (D) 
of figure 6. The initial choices of magnitudes ten and five 
for produced very slow convergence toward an admissible 
route, = 0.1, however, converged reasonably well. This 

can apparently be attributed to errors due to linearization 
on a path with relatively large curvature. 

The D-method produced rapid convergence to solutions for 
each of the last two cases. 

3.4 Conclusions 

On the basis of this study the following observations 
regarding the points of comparison (a) through (d) of 
section 3.1 seem justified. 
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(a) The choice of acceptable starting parameters is di- 
stinctly easier for the FL-method. The matrix p(t) is the 
set of actual control parameters for the problem and as such 
it is usually a set of well understood physical quantities,, 
These are less abstract than the adjoint vector or the con- 
stants Ch of the D-method. General knowledge about the 
physical concepts involved usually permits estimation of a 
starting control program which relatively closely approximates 
that of the optimum solution. 

(b) The interval of convergence, i.e., the domain of 

starting parameters from which a solution can be produced, is 
significantly larger for the FL-method. Any reasonable esti- 
mate for the control parameter will start the solution. Table 
I and figure 3 indicate there are seemingly good estimates 
for such as (-20,-80) which do not produce solu- 

tions for the D-method. 

(c) Once a workable set of starting parameters is found, 
however, both methods produce solutions from it for a wide 
range of variations of the problem. 

(d) If the D-method converges it consistently results in 
the more rapid convergence. For similar convergence criteria 
the computation time ratio is about one to two, and the number 
of iterations required, a ratio of about one to four. The 
criteria employed were different for the two methods, however, 
and it appears from table III that the requirement dg/ds < 10 
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was too tight for practical purposes, The precision is avail- 
able, however. 

The FL-method took only about ten percent more iterations 
to converge from a relatively wild estimate for p*Xt) than it 
took for a very close estimate. 

The following comments regarding the two methods are 
based on experience in implementing the two methods and on 
analyses of the papers by Faulkner [4] and Bryson and Denham 
[ 6 ]. 

The requirement of storing an entire control program and 
a complete set of functions for any one iteration is a 

disadvantage for the FL-method , This could result in memory 
storage problems for large problems or small computers. In 
contrast the D-method holds only the parameters C i and T or 
their equivalent. 

As indicated in Appendix B the FL-method has been develop- 
ed to the point of providing for variations at t=t Q and for 
the insertion of weighting functions by the matrix W(t), The 
power of these features has not been investigated in this 
paper. 

The problem of discontinuous control (corners) or un- 
bounded control have also been ignored. The D-method has been 
developed to handle them, [4], while the application of the 
FL-method to such problems has not been undertaken, nor is it 
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clear how this could be done. 

With the D-method there is no question of how closely we 
have approximated an extremal, with the mathematical routine, 
since only extremals are used. The only limitation is the 
accuracy of the integration method. Also, since only extremals 
are used various identities associated with an extremal can 
be employed. For example, if the independent variable does 
not enter explicitly, the Hamiltonian is constant. 

The integral relations for the corrections in the FL- 
method seem to be more stable from the computational standpoint, 
while, if the strong Legendre condition is not satisfied the 
routine must be modified in the D-method. 

Only forward integration (or only backward integration) 
is used in the D-method. In the FL-method the basic equations 
are integrated forward and the corrections are obtained by 
backward integration of the adjoint system. This necessitates 
either storing the complete trajectory just produced or re- 
constructing it from terminal values by backwards integration 
of the basic differential equations of the system. 
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APPENDIX (A) THE DIFFERENTIAL METHOD 



A.l Outline of the analysis 

The following relations have been defined: 



(3) 



x = f(x,p,t) the state variable equations. 



(4) 

( 12 ) 

(14) 



h( x, t ) ^ 



a = -F'A 



N terminal constraints, 



the adjoint system, 



r T 

[A'6x]m = / A' G 6p dt the fundamental 

Jt differential formula. 



( 21 ) 



x ( t ) = c 1 A 1 + + 



'2 '2 



c n A n 



n 



where A^(t) = 



An 

A 21 



A 



ni 



is the solution to (12) such that 






1, i = 4 

IP, i * J 



and x i = ° 

We write (14), using (21), as 

rT 

(A-l) [ X' 6x] t = X * ( t) G 6p dt 
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From (20) and (21) we get the m equations 
(A-2) A’(t) G(t) = 0 . 



It will be recalled that (12) with (A-2) constitute the 
Euler-Lagrange equations „ 

Additional relationships between the quantities dT, 
6x^, fip and dc^ will now be developed,, See [4], 

A differential dT in T results in the differential 

dx^(T) = 6x^(T) + x^dT , i = 1,2, « . . ,n , 



or 



dx T = fix^ + x dT 



giving the variations 



( A-3) 



fix = (dx - x dT), 



(A-3) with ( A-l) gives 

r T 

(A-4) X ' (dx - x dT) T = U'dx-Xx dT] T = / [ x'( t) G( t) 6p( t) ]dt, 

*0 

If x is held constant in (A-2), which in expanded form is 

I , .1 

X'G = c 2^2 + 000000+ c nr? = 0 > 

we can write 



(c l A ’]a 5 5p + A 'l° + 3 ° ° ° • ( c n A nfp 6p + A n° do n ) = 



0 , 
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where we define the m x n matrix 



8G, _ cKj , SCI , 

8p P ap-j_ P i p m 



m 



In terms of (21) this is the set of m 



equations 



(A-5) X §~6p + A' X G dc x + A' 2 G dc 2 + + A' m G dc n = 0, 



Solving (A-5)^for 6p^, 6p 2 , 6p m and substituting 

into (A-4) we get n equations, in the n+1 variables 
T,c 15 ,c n , of the form 



, n , 

( A-6) [ A^dx = 2 I ij dc j + dT 1 = 1 > * • • • » n 

3 ”1 



The differential of the constraint 



( A-7) 



dh k = 



^8^ dX i + 8t kdT ] T > 



For admissible curves h k is specified, 
therefore satisfy the relation 



k 1,2,.«.«,M • 
The end point must 



(A-8) dh^. = 0, k = 1,2, ,N. 

In vector notation (A-7) and (A-8) may be expressed by 

( A-9) [ Vh k » dr ] T = 0; k = 1,2,....,,N, 

where 7h k is a gradient in x,t space. 

For a maximum value of g(x,t),p we similarly have 

(A- 10) [ Vg»dr ] T = 0 

in x,t space indicating the differential of the terminal 
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point must be such as to result in displacements (geomet- 
rically speaking) normal to the Vg in x,t space, that is, 
tangent to the surface g(x,t) T = g ( maximum ) . 

Since the functions g(x,t) and h(x,t) are assumed 
to be independent the N+l by n+l matrix of coefficients 
from (A-9) and (A- 10) 
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9x n 
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has rank 


N + 


i. 










The 


transversal 


condition 


, [4] 


and 


be stated 


as 


the 


condition that the 


rank 



n+l matrix 



(A- 12) 



4 - 

M = 



X 1 X 2 



M 

• X 



n 



-X’x 



t = T 



be N+l and the rank of the matrices obtained by omitting 
either of the last two rows also be N+l, Row number N+2 
of M + is the row of coefficients of dx and dT of (A-4), 
These conditions on (A-12) yield n-N relations involving 
the terminal values of x, X, and t which must be satis- 
fied by a solution which gives at least a stationary value 
to g(x,t) T<! We designate them 
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s = 1.2 



n-No 



(A-13) H s (x,X,t) T = 0, 



^o*#**oooj 



The differential relations derived from (A-13), at the 
terminal point, are 



( A- 14) 




= 2 £! Is dx, + ZZ £2s£ii dC.+ fr^s jL + 

i Sx^ 1 i j ax i ac j j 'iaXj. 




dT 



T 



s = 1,2, . . . . . ,n-N, 



where X i is defined as in (21). 

Equations (A-6), (A-9) and (A-14) are 2n equations 
in the 2n+l differentials dx. , dT and dC., i and J= l,«,n. 



A. 2 Calculation procedure 

(a) A solution of the problem is characterized by the n+1 

parameters TjC^jC^, . . . 6 C n . We arbitrarily choose one 

v 2 

relation among the C*. We might choose iL(C.) = 1 or Just 

^ L 1 

designate one, say C n = 1. We must be careful only that the 
designated does not turn out to be negative or zero. We 
then guess a set of values for the remaining parameters and 
compute simultaneously a solution to the original differential 
equations (3) , the adjoint system (12) and the Euler 
equations (A-2) or (20). This yields an extremal E Q 
which does not, in general, satisfy (4) or (A-13). 

(b) Simultaneously calculate the correction integrals (A-6). 

(c) Treat the differential relations (A-5) through (A-13) 
as differences evaluated at the end point of E Q . We then have 
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( A“ 1 & ) A h^ ' h^ ™ h^ ^ ^ } ^0 J — 1 )»#«i**jN 

A H.^ "* H ^ ( -. ) , 1 j t * ( fl g fl e » I • jUl 

(d) Solve for the n differences AC-^, . . . . ,AC na>1 , AT and 
apply to the previous estimate for these parameters. If the 
original estimate was close enough the new parameters will 
result in a more nearly admissible extremal. 

(e) The process is repeated until the sequence of extremals 
converges to an admissible extremal within the limit of some 
convergence criterion. Additional analysis is required to 
distinguish between maxima, minima and stationary values. 
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Appendix B THE FL-METHGD 



Bel Outline of analysis 

We have previously defined the relationships 



(3) 


x = f ( x, p , t ) 


n state variable equations 


(4) 


h(x,t)^ — 0 


k constraints 


(10) 


6x = F ix + G 6p 


variational relationships 


(12) 


A= -f'A 


adjoint differential equations 


(14) 


T 

[/ / \fix] T = /\'(j 6p dt 


fundamental differential 
formula 



Following Bryson[6] we define one of the functions of 
(4) as a stopping condition, 



(B-l) -ft(x,t) T = 0 



which produces a terminal value T 
ing the need for guessing T« 

In section (1 0 5) we developed, 
adjoint system, 



(17) 



6h .( T) 

v 



-u 



tAh. G Sp dt > 



for t, thus eliminat- 



from solutions of the 



3 = 1 , 2 ,- 



,k 






or equivalently 



T 

(B-2) 6h(T) = f t /i h 



G 6p dt 
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where is defined as 



\<t>= 



A, 



h. 



• • • • 

H cn n to 



Au Ai_ 

n 12 n 22 



• • • • A i 



k2 



h 



In 



‘h. 



kn 



nxk . 



We also developed 



x 

(18) 6g(T) = /V G 6p dt . 

^ O © 



We now add, for (B-l) , a similar relationship 

T 



(B-3) 



M1(T) 






G 5p dt 



where A n is the solution to (12) such that 



A n (T) = V x A(T) , 



where 7 X denotes a gradient in x-space. 

For small perterbations the value of T will be changed 
only a small amount dT, so that 

dg = 6 g + gdT 
(b-4) dh = 6h + hdT 

da = 6 n + hdi 

.mere g, h, A, are defined as 
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* = (It + af f ) T 

W§* + &) t 

a = (at + a £ j i 



Accordingly (18), (B-2) and (B-3) modified for a change 
dT in T become 



(B-6) 



dg = 



dh = 



1 

Ag G 6p dt + g dT 
T 

f t A' h o 6p dt + h dT 



da = /\' a G 6p dt + n dT 



In order to limit the magnitude of the step 6p 
such that the linearization (10) is reasonable Bryson im- 
poses the limitation 

(B-7) (ds) 2 = 6p 6p dT 

o 

where (dS) is a magnitude arbitrarily chosen as a reasonable 

upper limit to the integral (B-7) „ 

Throughout this paper a simplified version of the FL- 

method, as developed by Bryson, has been employed in that 

variations at t = t have not been considered and in that 

o 

his matrix of 'weighting functions', W( t) , has been 
neglected. 
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ir equations (b-6) gives 



Eli’.: ina tie n c f c T 



(b-8) 


T 

dg = Jl AVn £ C *^ 

C 


(B-9) 


dh = A'A'ha 0 6p dt 




O 


v/here 





A sn = A g ' J- A ft A hfi = A h - ± h' . 

From (B-7) , (B-8) and (B-9) we form the linear 
combination 



T 

(B-10) dg = ( A' gn G -v’/^ G - U fip) 6p dt + v'dh +,/ (dS) £ 

The second variation of (B-10) is 

T 

«Ug) = X< A '?n 0 • vV 'hn 0 • £)j6p ' ) dt 

2 

For maximum dg the coefficient of 6 p = 0 implies 
(B-n) fP = </\ gq -A hfl v). 

Eliminating v and n between (B-7) , (B-9) and (B-ll) v/e 
have 



(B-12) 






6p(t)-± G 1 (/\ gn -A hn I hh 1 hgh J 



(d s) £ - dh'I'jj; dh 

_ T' T I 
■gg x hg 1 hh i hg 



+ G ' A ho dh 



-*■ Uso the plus sign for maxima. 
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where 




B.2 Calculation procedure for FL-method 

1. State a reasonable nominal control variable program, 
p*(t) which is defined for the time period (t c ,T). In 
many problems this can be a constant program. For ex- 
ample if the parameters are two directions, say the 
vertical angle <J) and the bearing 6, we could esti- 
mate p*(t)'= (<t> 0 ,© 0 ), where <j> 0 is any constant value 
such that 0 < <() 6 < tt/2 and © o is somewhere in the in- 
terval (correct heading ± tt/2) . Usually better estimates 
can be made. 

2. With p*(t) solve system (3) determining T and x(t) 
for t 0 < t v < T. 

3. Solve the adjoint system (12) backwards simulta- 
neously building A hn , A gfJ , 1^, I hg > I gg . Store 

^gn a ** ^hn G • 

4. On the basis of x(T) from step two, evaluate 
h(x,t) T o Set dh = -h(x,t) T ♦ 

5. Select a reasonable value for dS. If the numerator 
of the radical in (B-12) is negative scale down dh 
such that the numerator is ^.0. 
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6. Produce 6p*(t) from (B-12) for t t T 

simultaneously forming p (t) = p (t) + 6p (t) and 
store for the next iteration. 

7. Return to step two forming a new solution until the 
halt criteria are satisfied. 



Appendix C Flow diagram for adjustment of (ds) 

Definition of symbols 

dh -the error in meeting the constraints h(x,t)^ = 0. 

e-^ -the magnitude of dh within which it becomes 

reasonable to be concerned about extremality. 

Cg -the convergence criterion for extremality. 

-a part of (ds) whose magnitude is adjusted down as 

the solution approaches an extremal. 

D. -the denominator of the radical of (B-12) , 

1 1 

dg/ds in p- space. 




X D i < ^2 could also be used here. 
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